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We study the critical properties of the two-dimensional (2D) XY model in a transverse magnetic 
field with filling factors / = 1/3 and 2/5. To obtain a comparison with recent experiments, we 
investigate the effect of weak quenched bond disorder for / = 2/5. A finite-size scaling analysis 
of extensive Monte Carlo simulations strongly suggests that the critical exponents of the phase 
transition for f — 1/3 and for / = 2/5 with disorder are those of the pure 2D Ising model. The 
relevant low energy excitations are domain walls, and we show that their properties determine the 
nature of the phase transition. 
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I. INTRODUCTION 

In this paper we examine the frustrated XY model in 
two dimensions for two different values of the magnetic 
field representative of "commensurate states" . Exper- 
imental realizations of this model in the form of two- 
dimensional arrays of Josephson junctions and super- 
conducting wire networks [1-3] can and have been con- 
structed and one of the objectives of this work is to under- 
stand the results of these experiments. A perpendicular 
magnetic field induces a finite density of circulating su- 
percurrents, or vortices, within the array. The interplay 
of two length scales - the mean separation of vortices and 
the period of the underlying physical array - gives rise to 
a wide variety of interesting physical phenomena. Many 
of these effects show up as variations in the properties 
of the finite-temperature superconducting phase transi- 
tions at different fields. In recent experiments on super- 
conducting arrays the critical exponents of a number of 
these phase transitions have been measured [3], opening 
the opportunity to do careful comparison of theory and 
experiment. While we will discuss the model within the 
context of superconducting networks, the model is also 
closely related to the physics of adsorbed films on sub- 
strates which impose a periodic potential which differs 
from the preferred period of the adsorbed film. In this 
work we examine the ground state properties, low en- 
ergy excitations, and critical properties of the 2D XY 
model in the densely frustrated regime (/ 3> 0) for two 
particular values of the magnetic field. In addition, we 
investigate the effect of disorder on the ground state and 
critical properties. This paper elaborates and expands 
upon our previous results reported in Ref. [4] . 

The Hamiltonian of the frustrated XY model is 



u = - 22 % 



6j Ay), 



(1) 



where 6j is the phase on site j of a square L x L lattice 
and = (27r/</>o) A • dl is the integral of the vector 



potential from site i to site j with 4>q being the flux quan- 
tum. The directed sum of the Aj around an elementary 
plaquette ^2 Ay = 27r/ where /, measured in the units of 
0o i is the magnetic flux penetrating each plaquette due 
to the uniformly applied field. We focus here on the cases 
/ = p/q with p/q = 1/3 and 2/5. 

The ground state fluxoid pattern for these / is shown 
in Figure 1(a) [5,6]. The pattern consists of diagonal 
stripes composed of a single line of vortices for f = \ 
and a double line of vortices for / = |. These diagonal 
lines of vortices can sit on q sub-lattices and, in addition, 
there are q more states with the stripes going along the 
opposite diagonal for a total of 2q degenerate states. A 
common speculation for commensurate-incommensurate 
transitions and the frustrated XY model is that the tran- 
sition should be in the universality class of the q-state (or 
2q-state) Pott's model. We find that this is not the case 
because, as discussed below, domain walls between the 
different states vary considerably in both energetic and 
entropic factors. 

The effect of quenched impurities on phase transitions 
is an important and fascinating problem. The "Harris 
criterion" [7] indicates that the addition of (bond) ran- 
domness to systems which exhibit second-order transi- 
tions in the clean case with a positive specific-heat expo- 
nent a changes the numerical values of the critical expo- 
nents [8]. It has also been shown using phenomenologi- 
cal renormalization-group arguments that the addition of 
bond randomness to systems undergoing first-order tran- 
sitions results in a random-field mechanism at any coex- 
istence region which can cause the transition to become 
continuous [9]. Aizcnman and Wehr [10] have shown 
quite rigorously that in 2D a quenched random field re- 
sults, quite generally, in the elimination of discontinu- 
ities in the order parameter conjugate to the fluctuating 
field. Most cases where bond disorder has been studied 
and observed to change the order of the transition are 
for q-state Potts Models, where for q = 8 Chen et al. 
[11] found through extensive Monte Carlo simulations, 
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that the first-order transition of the pure model became 
second-order with the critical exponents being consistent 
with the universality class of the two-dimensional Ising 
model. Unlike q-state Pott's models with high q, the 
frustrated XY system is more readily compared to ex- 
periments such as recent experimental measurements of 
critical exponents in superconducting arrays [3]. 

II. STAIRCASE STATES 

The ground states of the Hamiltonian (1) will be 
among the solutions to the supercurrent conservation 
equations dTi/dQi = 0: 

^sin(^v-^-Ay,)=0 (2) 

y 

where f are the nearest neighbors to i. One set of so- 
lutions to these equations was found by Halsey [6] by 
considering the restriction to a quasi-one-dimensional 
case where one has adjoining staircases of current (see 
Fig 2(a)). All gauge invariant phase differences j m — 
On — @m — A m „, within a given staircase are equal and 
indexing the staircases by m as shown in Fig 2(a) one 
finds 

7m = 7r/m + a/2 — 7rnint[/m + a/(2ir)] (3) 

where nint is the nearest integer function, and a = for 
/ = p/q with q odd and a — ir/q for q even [6]. 

The staircase fluxoid pattern for / = 1/3 and 2/5 is 
shown in Fig 1(a) [5,6]. The pattern consists of diagonal 
stripes composed of a single line of vortices for / = ^ and 
a double line of vortices for / = |. (A vortex is a pla- 
quette with unit fluxoid occupation, i.e. the phase gains 
27T when going around the plaquette.) The stripes shown 
in Figure 1 can sit on q sub-lattices, which we associate 
with members of the Z q group. They can also go along 
either diagonal, and we associate these two options with 
members of the Z 2 group. In all, there is a total of 2q 
degenerate states(/ = p/q ^ 1/2). A common specula- 
tion for commensurate-incommensurate transitions and 
the frustrated XY model is that the transition should be 
in the universality class of the q-state (or 2q-state) Pott's 
model. We find that this is not the case because domain 
walls between the different states vary considerably in 
both energetic and entropic factors. 

III. DOMAIN WALLS 

Figure 1(b)- (e) shows the fluxoid pattern for some of 
the domain walls for / = 1/3. The domain walls can 
be classified into two types. Shift walls involve a shift of 
the vortex pattern across the wall (such as in Fig 1(b) 
where the pattern on the right is shifted down by one 
lattice spacings with respect to the pattern on the left) 



but the lines of vortices are still going along the same di- 
agonal. Herringbone walls are walls between states with 
the vortex stripes going along opposite diagonals. Note 
that there are q different walls of each type. 

These walls also have differing topologies. A herring- 
bone wall is very similar to a domain wall in an Ising 
model in that it separates two members of a Z2 group. It 
cannot branch into other herringbone walls and a 90 de- 
gree turn in the wall can be accomplished without chang- 
ing the vortex pattern, with the caveat that one considers 
the wall to be composed of sections of length equal to the 
distance between the diagonal lines of vortices (see Fig- 
ure 3). Thus, if one only has herringbone walls in the 
system, the set of possible domain wall configurations is 
similar to those in an Ising model. Shift walls, on the 
other hand can branch, both into other shift walls (with 
the constraint that the sum of the shifts on the walls af- 
ter the branch be equal to the original shift) and into a 
pair of herringbone walls, as shown in Figure 3(b). Shift 
walls also have an associated directionality in the sense 
that an attempt to make a 90 degree turn in a shift-by-n 
wall results in the wall changing to a shift-by- (q-n) wall 
(see Fig. 3(a) for an illustration). Since different shift 
walls can have quite different energies (see below) one 
finds that bends such as the one shown in Fig. 3(a) are 
energetically highly unfavorable as it can change a wall 
with low energy into a wall with a very high energy cost. 
A more energetically favorable kink in a shift wall can 
be formed by displacing a mismatched vortex on the wall 
in a direction parallel to the wall [12]. This displaces a 
section of the wall one unit cell in the direction perpen- 
dicular to the wall (see Figure 3(c)). Typically, one finds 
only kinks like these of size one or two lattice constants. 
Larger kinks start to produce long range distortions in 
the phase field and have higher energy. 

In order to calculate the energies of different structures, 
we solved the equations (2) numerically, using a quasi 
Newton method, on lattices with up to 2.3 x 10 5 sites with 
constraints fixing the fluxoid occupation of each plaque- 
tte (see Appendix A). Table I lists the energy per unit 
length a for straight domain walls between the various 
ground states at zero temperature for / = 1/3 and 2/5. 
One can see from the table that there is typically one or 
two walls with considerably lower energy than any of the 
others. Some of the patterns of energies seen in the table 
can be understood by counting the number of extra vor- 
tices in next or next-next-nearest neighbor plaquettes for 
the vortices along the wall. For instance, the energy of a 
/ = 1/3 shift-by-one wall is about twice that of the stan- 
dard herringbone wall. Looking at Figure 1(b)- (e) one 
can see that if you count the number of next-next-nearest 
neighbor vortices for vortices along each side of the wall, 
the shift-by-one wall has twice as many as the herring- 
bone. Similarly, walls which place vortices on nearest 
neighbor sites tend to be of a higher energy, or may not 
even be stable. While this does give a rough guide to the 
pattern of energies, it does not allow a strong comparison 
of walls with differently spaced vortices. 
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One can see from Figure 3(b) that a shift wall can be 
viewed as two adjacent, or bound herringbone walls. For 
/ = 1/3 the energy of two herringbone walls is less than 
that of a single shift wall and hence, the shift walls should 
be unstable to breaking up into herringbone walls. As a 
result, one expects that in the f — 1/3 case if the tem- 
perature is high enough for domain walls to enter the 
system, the herringbone walls should be the only walls 
present at large length scales. While for / = 1/3, her- 
ringbone walls are the only stable walls, this is not true 
for / = 2/5. For / = 2/5 it is energetically favorable for 
two herringbone walls to bind and form a shift-by-one or 
shift-by-three wall. This can lead to more complex do- 
main wall structures and has an important impact on the 
nature of the finite temperature phase transition. These 
issues will be addressed in more detail below. 

We also numerically calculated the energy of domain 
walls that are not straight. Figure 5 shows the energy of a 
square closed domains, formed from herringbone walls, of 
linear dimension L unit cells in a system of size 120x120 
with periodic boundary conditions. We see that to a very 
good approximation, the energy scales linearly in L. One 
can, however, work out some corrections to this linear 
dependence due to the change in the vortex density at 
the corner of the domains. For instance in Fig. 3(b) the 
vortices at opposite corners of the square domain have 
either an extra next-nearest neighbor vortex or a miss- 
ing next-nearest neighbor. From a distance, this gives 
a quadrapole moment to the domain. As the 3x3 do- 
main shown in Fig 3(b) is the basic building block of 
larger domains, one can conclude that larger differently 
shaped domains will not have a lower moment (i.e. they 
will be neutral and have no dipolc moment). The in- 
teraction of two such quadrapole domains at a distance 
x, large compared to it's size L goes like — 6 J e ff(L/x) 4 
(if one assumes an isotropic (which is not really true) 
interaction of two "corner" charges like J e f / In x) . In ad- 
dition, the self energy of a square quadrapole goes like 
-2 J eff ln\/2 + 2 J e ff InL. 

Figure 4 shows the interaction of some square domains. 
One sees that the quadrapole correction is measurable 
and fits the expected functional form quite well, but that 
the constants J e ff do not match what one would ex- 
pect from an isotropic calculation. In fact, the system is 
not really equivalent to an isotropic 2D Coulomb gas, in 
that the direction along the diagonal lines of vortices in 
the staircase state is not equivalent to the direction per- 
pendicular to the vortex lines. We have also calculated 
the energies of rectangular domains and some other less 
regular shapes and they have qualitatively similar (same 
functional form) behavior. 

The next question is whether or not the quadrapole 
interaction is likely to be relevant. One can use an argu- 
ment similar to that used to argue for a transition in the 
unfrustrated XY model. If you consider the interaction 
free energy contribution of the quadrapole interaction, 
it should contain the energetic part — A/r 4 and an en- 
tropic part — BT\n(-Kr 2 ) from confining the quadrapoles 



to have a separation less than r (One could do a more 
accurate calculation of the entropy but it will still have a 
lnr dependence). At the distances at which the —A/r 4 
form is valid, the lnr term wins all the time (at finite T) 
and hence one can argue that the quadrapole interactions 
should not be relevant. 



IV. SPIN WAVES 

At low enough temperature, domains should be small, 
and one is tempted to expand the energy about the 
ground state configuration. In this treatment, the pe- 
riodic character of the angles is neglected, but the exis- 
tence of long range order in the vortex lattice partly jus- 
tifies this method. The model is replaced by a so-called 
spin wave approximation which involves expanding the 
Hamiltonian to 2nd order in Ay , where 0y = Of- ' + An 

g(0) 



and 9\j is a ground state configuration: 



on 



(0) 



A 4 , 



ij kl 



d 2 n 

d6 kl d6i 



(0) 



(4) 



By definition, (dH/dOij)^ = and we just have a 
quadratic form. The free energy per site associated with 
(4) is 



T 



1 

1 



IniL 



| \ 'IX 



exp(-|E A 



d 2 n 



(0) 



1, (* JX ~ 1/2 
-— In det— 

(3 \ 2ir, 



(5) 



where J is the Jacobian matrix, J xx ' = d 2 H/d8 x d9 Dl 
The spin wave correlation function is 



G sw (xi,x 2 ) 

= (exp[i(A xl -A X2 )]) 



(6) 



x exp 



A X J X , X , A x , + i(A Xl - A X2 ) 



= exp ^-t^A(xi,x 2 ) t J 1 X(x 1 ,x 2 )j , 

where A(xi,X2) is a vector with +1 and —1 in positions 
xi and x 2 respectively and zeros everywhere else. 
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For the unfrustrated case [13] J is just the discretiza- 
tion of the Laplacian operator (i.e. 3<j) w V 2 0, where 
V 2 = d 2 (j)/dx 2 + d 2 <fi/dy 2 and the partial derivatives 
are replaced with a finite-difference formula d 2 <j)/dx 2 rts 
(<f>(x l+1 ,y i )~2cf)(xi,yi) + <p(xi- 1 ,y i ))/a 2 , and a is the lat- 
tice constant). As a result, ^X(x, x') T J _1 X(x, x') can 
be approximated by the Green's function for the Poisson 
equation, 

1 r 
9(r) = — In— , 
2tt r 

where r = o/(2v / 2e 7 ), a is the lattice spacing and 
7 = 0.577216 • • • is Euler's constant. This yields 

G st0 (xi,x 2 ) ~ exp(-l/27r/?ln(|xi - x 2 |/r ) 

= (r /|x! - x 2 |) W ( 7 ) 

So the correlation function of the spin wave fluctuations 
decreases according to a power law behavior. This alge- 
braic decay of the spin wave correlation function is bro- 
ken by the unbinding of vortex-antivortex pairs at the 
Kosterlitz-Thouless transition [13]. 

In the general frustrated case, J is not a discrete 
Laplacian. The question is, do we get something simi- 
lar? The 1 /r 4 interaction of the domains studied in the 
previous section suggests that we do. Figure 6 shows 
^A(x, x') T J~ 1 X(x, x') for / = 1/3 along a slice in the 
x-direction in a finite size system with periodic boundary 
conditions along the direction of the slice. The envelope 
of this curve is well described by the sum of two log- 
arithmic functions, lnx + ln(L — x) (where the second 
term comes from the periodic boundary conditions). In 
addition to this logarithmic part, there is a periodic os- 
cillation, coinciding with the underlying vortex lattice. 
In addition to this obvious oscillation, the phase of the 
oscillation depends on the initial x (The correlation func- 
tion is not just a function of (x — x')). The effect of this 
initial x dependent phase at long distances should not 
be important. However, distortions centered on nearby 
sites, and between rows of vortices can partially cancel 
due to this phase difference. There is also an anisotropy 
between the directions perpendicular and parallel to the 
diagonal lines of vortices. This anisotropy can, however, 
be removed in a continuum picture by rescaling the co- 
ordinates. 

This modified lattice "Green's" function also has an 
impact on vortex interactions. The presence of the log- 
arithmic part ensures that the overall flux balancing 
(/ = (rii) where rii is the vortex occupation of plaquc- 
tte i) is maintained. However the vortex interaction en- 
ergy should contain an oscillating component coinciding 
with the underlying vortex lattice. The effect of such a 
component is not entirely clear, especially as the ampli- 
tude of the oscillation does not decay away at large dis- 
tances. Conventional wisdom would suggest that as long 
as we still have the logarithmic interaction of vortices, 
they should still undergo a Kosterlitz-Thouless type of 



unbinding transition and an associated jump in the hc- 
licity modulus [13]. As we shall sec in the next section, it 
is not entirely clear whether or not this actually happens. 
It might be interesting to try to go through and derive 
the Kostcrlitz recursion relations with the oscillations as 
some sort of perturbation to see if it is relevant, although 
it seems unlikely to do anything but renormalize the core 
energies. 

V. / = 1/3 

The fluxoid pattern for the two lowest energy walls at / 
= | was shown in Figure 1(b) and (d). One can see from 
Figure 3(b) that a shift wall can be viewed as two adja- 
cent, or bound herringbone walls. For f = \ the energy 
of two herringbone walls is less than that of a single shift 
wall and hence, the shift walls are unstable and break 
up into herringbone walls. As a result, we confine our 
discussion of the f = \ case to the herringbone walls as 
other walls should not be present at large length scales. 
The energy cost for dividing anLxi lattice into two do- 
mains separated by a solid-on-solid (SOS) wall stretching 
from one side of the system to the other is 

H s ingie{z} = baL + ba^2\zk- Zk-i\. (8) 

k 

The height variables z k take on integer values (6 = 3 
is the shortest length segment). The partition func- 
tion, Z = X){ Zfe } exp(— 7i/T) can be evaluated either 
by the transfer matrix method or recursively (see Ap- 
pendix B) [14]. The interfacial free energy per column 
is T = Tln[e b,T / T tanh(6cr/(2T))]. The zero crossing of 
T gives an estimate of the critical temperature. Plug- 
ging in the values for the f = \ herringbone wall gives 
T c = 0.1 9 J, in remarkable agreement with the value 
T c = 0.22J found in the Monte Carlo simulations de- 
scribed below. 

Being similar to Ising walls, herringbone walls cannot 
branch into other herringbone walls, thus the set of pos- 
sible domain wall configurations is similar to those in 
an Ising model. We label the fraction of the system in 
state (s,j) as m s j, where s = ±1 denotes the member 
of Z 2 , and j = 1, 2, 3 denotes the member of Z 3 . Below 
the transition, one state (s, i) spans the system. On this 
state sit fluctuating domains, bounded by herringbone 
walls, of each of the states (— s, 1), (— s, 2), and (— s,3) 
in equal numbers; so the Z% symmetry is broken for the 
(s, j) states, but not for the (—s,j) states. As the transi- 
tion is approached from below, the domains occupied by 
the (— s, j) states grow, with smaller domains of the (s, j) 
states within them. At the transition, the Z 2 symmetry 
between the ±s states is restored and, as a result, the Z3 
symmetry for the (s, j) states is also restored. 

The Monte Carlo simulations used a heat bath algo- 
rithm with system sizes of 20 < L < 96. We computed 
between 10 7 and 3 x 10 7 Monte Carlo steps (complete 
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lattice updates) with most of the data taken close to T c . 
Data from different temperatures was combined and an- 
alyzed using histogram techniques [15] (see Appendix C). 

If the largest fraction of the system is in state (s,i), 
then we have three Ising order parameters, Mj = (m Si i — 
m- s j)/(m s j +m_ S; j), j = 1 • • -3. On average, these Mj 
are the same so we just take the average as M. To cal- 
culate the m a .i, we examine the Fourier transform of 
the vortex density pk± at the reciprocal lattice vectors 
k± = ±1) of the ground state vortex lattices. Start- 
ing from the definition of the Fourier transform, and us- 
ing the vortex states given above, one finds 

P — = m±i,i + m ±h2 e t2 ^ 3 + m±i, 3 e- i27r / 3 , (9) 
Pa 

where p g is the modulus in the ground state. In practice, 
Pk± is reduced by small short-lived regions which don't 
quite match any of the six states. Since this effect is the 
same for all states, it cancels when calculating M. Us- 
ing the real and imaginary parts of pk± in addition to 
"Ylij m ±i,j, calculated from the direct vortex lattice as in 
[16], we can find the five independent m a j. 

In addition to the energy and order parameter, several 
other quantities were calculated from the Monte Carlo 
data using the corresponding fluctuation-dissipation re- 
lations: 



X = KL 2 ((M 2 )-(M) 2 ), 
«91n(M") _ (M n E) 



dK 



(M n ) 



(10) 



where K = J/ksT. In addition to the discrete order 
parameter, we also followed the helicity modulus defined 
by Y x , y — d 2 J r /d(j) 2 \ < p — 0, where T is the free energy 
density and <fi is a twist in the boundary condition along 
the x or y direction. The helicity modulus also follows a 
fluctuation-dissipation relation which is used in calculat- 
ing it from the data: 



Y * = J2\ E[( r - r ')-x] 2 COs(fl r 



\{r,r<) 




[(r - r') • x] sm{6 r - ( 

(r,r'> 

E I( r - r ') • *] sin(e r - e r 

(r,r'> 



where (r, r') denotes nearest neighbor pairs. 

To determine the critical exponents for the transition 
we make use of finite size scaling [17] . Following standard 
arguments, one assumes that for a second-order transi- 
tion, the singular part of the free energy, F(t, h), near the 




transition is dominated by a term that changes under a 
change of scale according to the ansatz 

F(t,h) = XF(X s t,X r h) 

where t = (T — T c )/T c and h is an applied field which 
couples to the order parameter M (so h is not the true 
magnetic field here). From this, one can derive the scal- 
ing form of the order parameter, specific heat, suscepti- 
bility, etc. using the standard relations, M — —dF/dh, 
C = -Td 2 F/dt 2 , X = dm/dh, etc. If one takes the 
special case h = 0, A = |t| -1 / s one can relate r,s to the 
standard exponents a for the specific heat, (3 for the order 
parameter, and 7 for the susceptibility as s = l/(a — 2), 
r = (7 + (i)/{a - 2) and a + 2/3 + 7 = 2. If one takes the 
case h = and A = £("~ 2 )/ l/ \ where v is the exponent 
for the divergence of the correlation length, one obtains 
the relations for finite size scaling: 

M n = L-M"M n (x t ), 
C = L a ^C(x t ), 

X = L^X(x t ). (12) 

where x t — tL x l v is the temperature scaling variable. 
Using relations 12 one can also derive [11,18] 



d(M) 
dT 
dln(M) 
dT 



= L—V(x t ), 
= L 1 ^Q(x t ). 



(13) 



For a finite lattice the peak in, for example the specific 
heat, scales with system size like C max oc L a / U and oc- 
curs at the temperature where the scaling function C(x t ) 
is maximum so that 



C(x t ) t 
dx t 



\x t =x* 



= 0. 



This defines the finite-lattice transition temperature 
T C (L) by the condition x t = x^ so that T C (L) = 
T c + T c x\L~ x l v . In general the finite-lattice transition 
temperature calculated from different quantities differs 
slightly but extrapolates to the same T c in the limit of 
large L. 

A very accurate way of locating the transition temper- 
ature is by using Binder's cumulant [19], 

U = 1- (M 4 )/(3(M 2 ) 2 ), 

shown in Figure 7. For system sizes large enough to obey 
finite-size scaling, this quantity is size independent at the 
critical point. From Fig. 7 we find T c = 0.2185(6) J. T c 
can also be determined from the scaling equation for the 
temperature at the peak of thermodynamic derivatives 
such as the susceptibility, T C (L) =T C + aL^ 1 ^. We find 
these other methods give T c in agreement with that from 
U. 
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Finite size scaling [17] at T c applied to d In M / dK gives 
\jv = 1.011 ± 0.029, and to the susceptibility \ gives 
7/i/ = 1.758 ± 0.013, and to M gives 0/u = 0.14 ± 0.02. 
These exponents are determined from the slopes of the 
lines shown in Fig. 9 which plots the values of these quan- 
tities at the critical point as function of L. These expo- 
nents are in excellent agreement with the Ising values 
v = 1, 7 = j, and = |. Fig. 8 shows the collapse of 
the raw data onto the scaling function (inset) for \- 

Two previous examinations of the f = \ case [12,20] 
suggested a continuous transition but did not measure 
critical exponents. Lee and Lee [16] claim to find sepa- 
rate, closely spaced transitions, for the breaking of Z 2 and 
Zz- One explanation for their conflicting results comes 
from the small system sizes (L < 42) used in their anal- 
ysis. Below the transition, if the dominant state is (s, i), 
in small systems you often do not see all three of the 
(— s, j) states in the system at the same time. Figure 
10 illustrates this effect. The minimum of (pk+,Pk-) is 
a measure of the Z% symmetry breaking for the (— s,j) 
states and this goes to zero as L — ► oo. The finite value 
of mm(p k+ , Pk-) for small L can give the impression of 
separate transitions for small systems (If a measured pa- 
rameter contains a contribution from vaim(p k+1 p k -) it's 
derivatives can have a double peaked structure from the 
derivative of mm(p k+ , Pk-))- One must take care in the 
choice of order parameter to ensure that this contribu- 
tion is not biasing the results. For example we found 
that the derivative of the Ising order parameter used in 
[16], M' — mi + m 2 + m 3 — m 4 — m 5 — m 6 has a double 
peaked structure for intermediate lattice sizes that does 
not completely go away until L = 96. This makes M' an 
unsuitable choice of order parameter for finite-size scal- 
ing. This is also the probable cause of the presence of a 
shoulder in the specific heat at intermediate system sizes 
[16]. For larger L, we see this shoulder merge with the 
main peak and for L = 84 and 96 it is no longer clearly 
discernible (see Fig. 11). 

The helicity modulus Y is the quantity most closely re- 
lated to experimental measurements [13]. For / ^ 0, the 
scaling of the I-V curves found in experiments is consis- 
tent with domain wall activation processes [3] . The the- 
ory of Nelson and Kosterlitz for the / = case predicts 
that Y should come down in a characteristic square-root 
cusp and then jump with a universal value, 2kgTxT/ 7T - 
However, we find an exceptionally good fit (Fig 12) of our 
data to Y-Yq = L-P/ v M((T -T C )L V I V ) with v = 1, 
— |, and Yo — 0, which is the scaling form of M. Clearly, 

Y is affected strongly by fluctuations in M and attempt- 
ing to fit scaling relations for the / = case [16] without 
taking this into account seems questionable. We see two 
possible interpretations of our result. The first is that 

Y only receives contributions from the ordered part of 
the lattice. So comparisons with the / = case should 
examine Y m — Y/M. Y m w 0.58 at the transition im- 
plying a larger than universal jump. Alternatively, one 
can say that although Y is brought down by fluctuations 
in M, it should still jump when it crosses the universal 



value, IkuT 1^- Extrapolating the observed behavior of 
Y gives Fl^oo = a\T — Tc] 13 . This crosses the value of 
the universal jump at Tkt — T c w 10~ 6 . Although we 
do not see evidence for a jump, a difference in transition 
temperatures of 10~ 6 would not lead to any observable 
effects for the system sizes studied here. 

VI. / = 2/5 

While for / = i , herringbone walls are the only stable 
walls, this is not true for / — |. For / = | it is ener- 
getically favorable for two herringbone walls to bind and 
form a shift-by-one or shift-by-three wall. Binding does, 
however, have an entropic cost. To see if these walls are 
bound we consider the following model for two SOS walls: 

H d {A, z} = ^2{(2b<r + U||^,o) + b<r\z k - z k -i\ 

k 

+ (2ba + u ± 5 Zk ,o)&k + V r ({A, z})}. (14) 

z k is the separation of the walls (z k > 0), A k is the 
number of vertical steps the two walls take in the same 
direction in the k'th column (— oo < < oo). uh and 
uj_ are the binding energies parallel and perpendicular to 
the wall. At this stage we take V r = 0. The solution to 
such a model is discussed in Appendix B. A ground state 
eigenvector ipn( z ) — where 1/p is the localization 

length, or typical distance separating the lines, charac- 
terizes the bound state of the two lines, p = defines 
the unbinding transition at T^. For the cases of inter- 
est, one finds = 0.398J for the shift-by-one walls and 
Xf, = 0.442 J for the shift-by-three walls. In addition, the 
free energy for these walls crosses zero before they un- 
bind. Hence, at the transition, defined by the point at 
which the walls enter the system, we expect a branching 
domain wall structure similar to the q > 5 Pott's models 
where a first order phase transition occurs. Technically, 
this is a mean field argument for the interfaces but, since 
the interfaces are extended objects it should give a rea- 
sonable picture of the order of T\, for the interfaces and 
T c . 

In their Monte Carlo simulations, Li and Tcitcl [21] 
observed hysteresis of the internal energy when the tem- 
perature was cycled around the transition and used this 
as an argument for a first order transition at / = |. 
The most direct indication of a first order transition is 
the presence of a free energy barrier between the ordered 
and disordered states which diverges as the system size 
increases [22]. The free energy as a function of energy 
is obtained using Tl (E) = — In P L (E) where P L (E) is 
the probability distribution for the energy generated by 
Monte Carlo simulation of a L x L system. Figure 13 
shows the growth in this barrier as the system size in- 
creases from L = 20 to 80 giving clear evidence for the 
first order nature of the transition. 

Since there is no diverging characteristic length to 
which the linear dimension L could be compared at a 
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first order transition, one finds that it is simply the vol- 
ume L d that controls the size effects [23] . One thus finds 

Cmax: X ^ L 

for a first-order transition. Figure 15 shows the spe- 
cific heat as a function of L 2 for the / = 2/5 clean 
system. The linear fit (solid line) clearly shows the ex- 
pected first-order scaling behavior. Similar behavior can 
be seen in the susceptibility as shown in the Figure. From 
the positions of the peaks as a function of L we obtain 
T c = 0.2127(2) J. 

VII. DISORDER AND THE / = 2/5 PHASE 
TRANSITION 

We now consider the effects of disorder on the / — | 
phase transition. Taking the couplings in the Hamilto- 
nian (1) as Jy = J(l + ey), the are chosen randomly 
from a Gaussian distribution with a standard deviation 
S. Due to variations of the phase differences across the 
bonds, a specific realization of random bonds may favor 
a certain sub-lattice for the ground state, creating an ef- 
fective random field. To quantify the effect, we placed 
the fluxoid configuration of the ground states down on 
10 000 separate realizations of the disorder and allowed 
the continuous degrees of freedom (the phases) to relax 
and minimize the energy. We find that the changes in 
energy from the 5 = case fit a Gaussian distribution 
with mean — 0.5S 2 L 2 and standard deviation 5L. The 
difference in energy between states which were degener- 
ate in the clean system is the measure of the random 
field. This difference centers on zero and has a standard 
deviation of 0.75<5L for two states related by a shift and 
0.575L for two states with vortex rows along opposite di- 
agonals. The effect of random fields on discrete degrees 
of freedom in 2D is marginal [24] . For D > 2 there is a 
critical randomness above which random fields cause the 
formation of domains in the ground state of size <~ £ r f. 
Aizenman and Wehr have shown that this critical ran- 
domness is zero in 2D [10]. Yet, their result does not 
preclude the possibility that £ r f is so large as to be unob- 
servable in a finite sized sample. Indeed, experiments on 
superconducting arrays have found apparent phase tran- 
sitions, including scaling behavior [3] in sample sizes of 
order 1000 x 1000. In our simulations with disorder at 
S < 0.1, all systems had a low temperature state with the 
order parameter approaching unity. We will, therefore, 
ignore the effects of random fields for 6 < 0.1 assuming 
that £ r f is larger than the sample size. 

At any coexistence point of the clean system, random 
bonds result in different regions of the system experienc- 
ing average couplings slightly above or below the critical 
coupling. As a result, at any given temperature the sys- 
tem will predominantly prefer cither the ordered or dis- 
ordered state wiping out the coexistence region and leav- 
ing only a continuous transition [24,9,10]. It has been 



conjectured [11] that critical random Potts models are 
equivalent to Ising models. Kardar et al. [25] suggested a 
possible mechanism for this effect. Their position space 
renormalization group approximation suggests that the 
probability of loop formation in the fractal interface of 
the clean system vanishes marginally at a transition dom- 
inated by random bonds. The interface may have some 
finite width due to a froth of bubbles of different phases, 
but under renormalization a linear critical interface is 
obtained and, hence, an Ising transition appears. 

The fluxoid configurations from our simulations sug- 
gest that for large enough disorder, (<5 > Sf) the interface 
is really linear, not just in the renormalized sense. Sf can 
be estimated by placing a random potential V r in Eq. 14. 
Ignoring the terms involving one obtains the model 
for wetting in the presence of disorder, solved by Kardar 
[26] in the continuum limit. He obtained a new length 
scale due to randomness, 

1/k = 2T 3 /KS 2 

where K is the renormalized stiffness related to the inter- 
facial free energy <j{6) by K = a(0) + d 2 a(8)/d6 2 \o where 
9 is a small tilt angle of the interface. For an Ising-likc in- 
terface K w T/asinh[kr/T - lncoth(&er/(2T))] [14]. The 
unbinding transition is lowered and is now defined by the 
condition [i — n = 0. As decreases, it eventually hits 
the transition temperature for the first order phase tran- 
sition observed in the clean system. At this point any 
branched domain wall structure is unstable. This is just 
the last step in a process in which the effective linear in- 
terface becomes narrower as disorder increases. In the 
vicinity of this "final" (mean-field) unbinding, the Ising- 
type behavior of the system should be readily visible at 
any length scale. 

We have done a Monte Carlo analysis with bond dis- 
order values of 5 = 0.05 and 0.1. Since we are dealing 
with quenched disorder, we are interested in averaged 
quantities; for instance the free energy is 

T= -k B T[\nZ} av (15) 

where the square brackets indicate an average over dif- 
ferent realizations of disorder. Since most quantities of 
interest involve derivatives of the free energy, to calculate 
the average value of a thermodynamic quantity, we first 
calculate it for a given realization of the disorder and then 
do a configurational average over 10 to 15 realizations for 
S = 0.1 and seven realizations for S = 0.05. Figure 14 
shows the free energy barrier for / = ■jr clS cL function of 
system size in the for S = 0.05, and 0.1. For S = 0.05, the 
barrier first grows with system size and then levels off. 
At S = 0.1 the free energy barriers are essentially zero, 
indicating a continuous transition and that the system 
sizes are large enough to apply finite size scaling. Here, 
we follow the finite-size scaling methods used in [11]. 

Figure 16 shows the peak values of d\nM/dK and 
X as a function of L. The slopes of these plots give 
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1/u = 1.05(12) and 7/1/ = 1.70(12). A similar analy- 
sis of dM/dK gives (1 - 0)/u = 0.94(10) [4]. As in the 
/ = 1/3 case, the helicity modulus appears to track the 
order parameter M. Within errors, these exponents are 
what one would expect from an Ising model. Experi- 
ments at / = I [3] also found a continuous transition 
and measured the critical exponents v = 0.9(5) and the 
dynamic critical exponent z = 2.0(5), consistent with an 
Ising transition. 

VIII. CONCLUSIONS 

In conclusion, we find that the nature and universality 
class of the phase transitions are quite sensitive to the 
proximity of the binding transition for the lowest energy 
domain walls. For / = 1/3 the lowest energy walls are 
never bound and the transition is Ising- like. For / = 2/5 
domain walls can lower their free energy by binding to 
each other, resulting in a first order phase transition. Dis- 
order weakens this binding and changes the transition to 
be continuous and Ising-like. These results are consistent 
with the continuous phase transition and critical expo- 
nents observed experimentally for f — 2/5 [3]. 

We thank M. Aizenman, P. Chandra, J.M. Kosterlitz, 
X.S. Ling, and D. Huse and for useful discussions. 

APPENDIX A: CONSTRAINED OPTIMIZATION 
FOR VORTEX LATTICES 

Minima of the Hamiltonian (1) satisfy Equations (2). 
However, these equations are written in terms of the 6j 
variables and the locations of the vortices does not enter 
explicitly. This is quite inconvenient as one finds that the 
zero temperature energies of the system are almost en- 
tirely dictated by the vortex structure. By this we mean 
that given the position of all the vortices, the phases 
appear to be uniquely determined (up to an overall con- 
stant) by the minimization conditions. This can be made 
more explicit by working with the gauge invariant phase 
differences 

2tt f&fi 

7»,j = - - — / A • efl, 

H>o J(i,j-i) 

2tt 

ai,j = 0i- hj - e id - — I A • dl, (Al) 

00 J{i,j) 

where 6ij is the phase on the site at row i column j of 
the lattice. This introduces an extra variable per site 
(instead of just 9ij now we have 7^ and caj) and a 
compensating constraint that 

lij - li-ij + a i,j - a i,j-i - M/ - riij) = 0. (A2) 

That is to say, the sum of the gauge invariant phase dif- 
ferences around any plaqucttc must equal the magnetic 



flux through the plaquette 27r/, plus an integer multi- 
ple riij of 2ir. If the gauge invariant phase differences 
are restricted to a range of 2tt such as [—n, n) then riij 
measures the vortex occupancy of the plaquette and is 
typically or ±1 with the sign depending on the sign of 
/• 

One then rewrites Equations (2) in terms of the gauge 
invariant phase differences to get 

sin 7^ — sin7j ; j + i + sinaj+ij — sinajj = 0. (A3) 

If disorder is added, the random couplings should be in- 
cluded here. These, in addition to Eq.'s (A2) give 2MN 
equations (for a lattice of M x N unit cells) for the 2MN 
unknown gauge invariant phase differences. The vortex 
pattern {riij} is now an input and stays fixed. When 
periodic boundary conditions are imposed one finds that 
two of these equations are not independent. Two more 
convenient conditions to impose closure are 

M 

sin a N j - I c = 0, 

N 

5^sin7i,i-J r = 0, (A4) 

i=l 

where I c is the net current flowing down the columns of 
the lattice and I r is the net current flowing along the 
rows. In all cases found, the lowest energy state corre- 
sponded to I riC = 0. 

The above equations can now be organized into the 
form F({jij, ctij}) — as 

AT 

Fi = ^shi7 4! i - I r , 

i=l 

F 2 M(i-l)+2j-l = - 

+a id - a i; j_i - 2n(f - n itj ), 
F 2 M(i-i)+2j = sin 7; j - sin7 ij+ i 

+ sina i+ ij - sma id , 

M 

F 2 mn = sin aN <i ~ Ic ' ^ A5 ) 
i=i 

If we define x to have elements x 2 M(i-i)+2j-i = 7»,j an( i 
x 2M(i-i)+2j = a i.j (i = 1- ■ ■ N and j = 1 • • • M) then the 
solution to (A5) can be found using Newton's method 
which involves iteratively solving 

J • <5x = F (A6) 

and updating x, 

Knew = *old + Sx, (A7) 

where the Jacobian Jij — dFi/dxj. 

The set of equations (A6) can be very large (we solved 
systems with up to 2.3 x 10 5 sites which means Eq.(A6) 
represents about half a million simultaneous equations) . 
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In addition, we need to solve these systems very fast, es- 
pecially when disorder is added and averages over tens of 
thousands of solutions are needed. This is made possible 
by the special form of the Jacobian matrix: 



J = 



(A8) 

where the dots represent the non-zero elements. We see 
that J is very nearly band diagonal. In fact J can be 
written as 



J = A + U • V T 



(A9) 



where A is the band diagonal part of J (the same three 
matrix diagonal blocks as J) and U and V are N x 2M 
matrices (as opposed to 2MN x 2MN). I should point 
out here that the method described below has a speed 
that is proportional to NM 2 so that the axes of the lat- 
tice should always be chosen so that M < N for efficient 
operation. U and V have the form 



(A10) 

The first two blocks of U and V T have the nonzero ele- 
ments indicated and and the remaining blocks of U are 
from the first block column of J and the remaining blocks 
of V are from the first block row of J. 

The solution of a band diagonal system A • x = b is 
considerably simpler than solving a general linear system 
of 2MN equations. Not only that, but the LU factoriza- 
tion of A has the same storage requirements as A which 
can be stored in a packed storage scheme holding only the 
central nonzero band. In order to solve our slightly more 
general problem we make use of the Woodbury formula 
[27]: 

J- 1 = (A + U- V T )- 1 
= A- 1 



- [A" 1 • U • (1 + V T • A" 1 • ITT 1 • V T • A" 1 ] . 

(AH) 

Since storage of A -1 is not practical (the inverse does 
not preserve the band structure of the matrix), we must 
make use of (All) in the following way, as described in 
[27]: To solve the linear equation 



(A + U • V T ) • <5x = -F 
first solve the 2M + 1 auxiliary problems 
A Z = U, 

and 

A • y = — F. 



(A12) 



(A13) 



(A14) 



This can be done by LU factorizing A once and then 
using the factorization to solve all the systems simulta- 
neously. Routines from LAPACK [28] can make this very 
fast and efficient. Next, do the 2M x 2M matrix inversion 



H= (1+V 1 



(A15) 



In terms of these quantities, the solution is given by 

<5x = y - Z • [H • (V T • y)] . (A16) 

In order to start Newton's method, one needs a good 
initial guess. This is provided by patching together the 
staircase state solutions described in section II. In ad- 
dition, care must be taken to ensure that the gauge in- 
variant phase differences do not wander out of [—n,ir). 
There are a number of options one can use if a phase 
difference wanders out of range. One is to just pin the 
solution at ±7r. This is not a great solution as this is 
not really a minima of the unconstrained Hamiltonian. 
Another solution is to just add or subtract 2ir and con- 
tinue iterating Newton's method. This can cause a jump 
in the errors on one of the equations which may result 
in a large change in x at the next step which may or 
may not be beneficial. Another solution is to replace the 
phase difference with the value on the other branch of 
the arcsin function on [— n, it). This causes no change 
in the error on the current conservation equations and 
produces a smaller change in the corresponding Eq.(A2). 
Many of these problems can often be avoided by taking a 
step in the Newton direction but with smaller length, es- 
pecially in the initial stages, using a dynamic step length 
algorithm similar to those described in [27]. 



APPENDIX B: SOLID ON SOLID MODELS 

A good review of interface models is given in [14]. Here 
we briefly discuss the cases relevant to our situation. The 
SOS model of an interface ignores overhangs and bubbles 
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and configurations can be described in terms of integer- 
valued height variables whose values are measured from 
the T = position of the interface (see Figure 17). The 
energy cost for dividing anlxl lattice into two domains 
separated by a solid-on-solid (SOS) wall stretching from 
one side of the system to the other is 



H s ingle{z} = bcrL + b<T ^ \ z k ~ Z k - 



(Bl) 



The height variables z k take on integer values (b is the 
shortest length segment). The partition function, Z = 
J2{ Zk } ex P( — 7~L/T) can be easily evaluated by change of 
variables, Aj = — so that 



Z=l[e ^ £ 



-/3b<jA k 



k=l 



A fc =- 



where [— r, r] is the allowed values of A k . In the unre- 
stricted case r = oo, the interfacial free energy per col- 
umn is T = Tln[e /3i " T tanh(6cr/(2r))]. The zero crossing 
of T gives an estimate of the critical temperature. In the 
case of the two-dimensional Ising model this zero crossing 
gives the exact critical temperature. This is somewhat 
fortuitous, but nevertheless useful. 

In the continuum limit, the problem of two interfaces 
can usually be broken down into a center of mass part 
and an independent part involving the separation of the 
two interfaces. We would prefer, however, to work with 
a discrete model with parameters input from the energy 
calculations of the appropriate bent domain walls. We 
were unable to find the solution to such a model in the 
literature, so we present one here. Questions that we are 
interested in are whether or not the two interfaces are 
bound and whether or not unbinding occurs before or 
after the free energy of the walls becomes negative. To 
answer these questions wc consider the following model 
for two SOS walls shown in Figure 18: 

Hdoubie{A, z] = ^{(2kr + U||(5 Zfci0 ) + ba\z k - z k -i\ 



+{2ba + u ± S Zk<0 )A k }. 



(B2) 



where z k is the separation of the walls (z k > 0), and A k 
is the number of vertical steps the two walls take in the 
same direction in the k'th column (— oo < A k < oo). uh 
and u± are the binding energies parallel and perpendic- 
ular to the wall. 

The partition function is 

L 

Z =^2Y[ e~ l3ba \ Zk ~ Zk - 1 \e~ l3( - 2 '" J+u u Sz f ' > 

{z k }k=l 

x{(l + \z k -Z k - 1 \) + ]T e -f^ + u ± S^ )lA k \ } 



A fc #0 



(B3) 



The (1 + \z k — Zk-i |) comes from the fact that for A^ = 
there are \z k — Zk-i\ + 1 ways to divide the change 



\z k — z k —i\ between the two lines. Summing over A k 
leaves the partition function in the form of a transfer 
matrix: 



L 



-fiba\z k -z k - 1 \ 



{$z k ,o{\zk - Z k -l\ 



{z k } k=l 



+ coth[(3(ba + u ± /2)})e-^ 2ba+u ^ 

+ (1 - S Zkfi )(\z k - « fc _i| + coth/3He- 26,T } 



{z k }k=l 



(B4) 



Unfortunately, we were unable to solve the general case 
analytically. However, restricting z k — z k -i to or ±1, we 
can derive the eigenvalues and eigenvectors of the matrix 
T explicitly. A ground state eigenvector ip^{z) = e^^ z , 
where 1/fj, is the localization length, or typical distance 
separating the lines, characterizes the bound state of the 
two lines. Vv( z ) is found by first finding the eigenvalue 
Xfj_ (from the defining equation (Tip^z = X^ip^z)) for 
z > 0. ji is then obtained from the eigenvalue equation 
for z = 0. This gives e M as the solution to the quadratic 
equation, 

(1 + coth/%cr)e 2 ^ 



+e 



f3ba 



coth [3ba - e- pu w (1 + 2e-^ 2brT+u ^) 



1 + coth [3ba - 2e- /3u " (1 + e -0Wv+u ± )-j\ = q ( B5 ) 



ji — defines the unbinding transition at Ty,. The more 
general case, \z k — z k -\\ < N with N a large number (typ- 
ically about 1000), can be easily solved numerically and 
is not that different from the restricted case discussed 
above. The values quoted in the text arc from such a 
numerical calculation. 



APPENDIX C: MONTE CARLO SIMULATION 
OF CONTINUOUS SPIN SYSTEMS 

A reasonable introduction to Monte Carlo techniques 
is given in [29]. However, some of the implementation 
techniques suggested in this book are out of date and 
should be taken with a lump of salt. Most simulations 
of frustrated spin systems described in the literature ap- 
pear to have used a rather poor updating scheme leading 
to very long autocorrelation times. We use a heat bath 
scheme described below which seems to be a couple of 
order of magnitude faster than these standard schemes 
near the critical point. This is not to say that other heat 
bath schemes have not been used, it is just that such 
works almost never describe any details of how this is 
done, a problem we shall try to rectify here. To make 
efficient use of the data generated in a Monte Carlo sim- 
ulation one should make use of the histogram techniques 
of References [15,22]. 
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1. Sampling 

Formally, the task of statistical mechanics is to com- 
pute from the model Hamiltonian Ti the desired average 
properties, 

<^({%})> = \\ 4%M«M)cx P [-w({M)/t], 

(Cl) 

where states are weighted with the normalized Boltz- 
mann distribution 

P({%}) = ^exp[-W({%})/T]. (C2) 

While this gives a formally exact description of the 
probability distribution, we are not really interested in 
such detailed information, nor is it possible to carry out 
the integrations in the high-dimensional space required 
in the thermodynamic limit. The dimension of the space 
can be reduced somewhat by making use of finite size 
scaling to extrapolate from small systems (L < 100) to 
the thermodynamic limit. Even for these smaller L, it 
is still not possible to numerically integrate the system 
based on any sort of discretization scheme. One instead 
uses Monte Carlo integration which is simply to pick N 
sets of {Oij} randomly distributed according to (C2) and 
then 



1 N 



(C3) 



1=1 



If the {&ij}i are independent and A({0ij}) is distributed 
in a Gaussian distribution with variance a 2 then the error 
in (A) calculated in this manner is er/A 1 / 2 . 

In practice, the knowledge of how to pick independent 
random numbers distributed according to (C2) is quite 
close to knowing how to solve the problem exactly. In 
general, we must give up on the idea of independent 
random numbers and instead construct a Markov pro- 
cess where each state {9ij}i+i is constructed from a pre- 
vious state {%}; via a suitable transition probability 
W({6ij}i — ► {8ij}i + i). A sufficient condition for the 
distribution function P ({Oij}) of states generated to con- 
verge to (C2) in the limit N — > oo, is for the transition 
probability to satisfy detailed balance: 



wMAv — {%}*) 



= exp - 



sn 



(C4) 



where SH = — Note that equation 

(C4) must be satisfied for all possible moves / — > I' in 
order to be ergodic. 

This still leaves many choices for the move. Ideally, 
one would like to change many degree's of freedom si- 
multaneously, unfortunately in the absence of any clus- 
ter routines for frustrated systems, one is left with single 



site updating moves. (Alternatively one can simulate a 
Langevin equation to change all degree's of freedom si- 
multaneously, but by a small amount. Even Langevin 
dynamics are not unique, and the dynamics which are 
supposed to be appropriate for superconducting arrays 
[12] was found to have longer autocorrelation times than 
the Monte Carlo method we ended up using.) One partic- 
ularly poor, but popular, method of updating continuous 
degrees of freedom involves picking a new Oij completely 
at random, or in an interval about it's previous value, and 
then accepting or rejecting the move based on whether 
another random number is above or below exp (— 
This can give extremely long autocorrelation times, and 
leads to a high number of rejected moves in the low tem- 
perature state. One would have to apply this same step 
numerous times to the same spin just to equilibrate it 
with it's nearest neighbors. 

An ideal single site updating step would pick Oij ac- 
cording to the conditional Boltzmann probability p(0ij) 
for Oij given the knowledge of the neighboring spins 
{0ij±i,0 i± i j}. For our frustrated XY model this is 



,.7 + 1-1 



= £ CX P [( cos ( e i,j+i - °i] + A ij 

+ cos(6> i;j _i -0^ +A*/ _1 ) 
+ cos(0ij-0 i+ld +AV +ld ) 
+ cos (0ij-0i- 1 ,j+At hj ))IT 

= j^jhj CX P (jp cos (% - 5 )) . 



(C5) 



where 



h = \J x 1 +y 2 , 
S = arctan(x/y), 

x = MOu+i + + M0ij-i + 4- _1 ) 

+ - AV +l j) + sin(^_ 1;j - AV_ hj ), 

y = cos(0 <J+ i + A^ +1 ) + cos(^-_i + A^- 1 ) 

+ coa(e i+ltj - A% hj ) + cos^-! j - A£ ld ), (C6) 

and Iq(x) is the zeroth order modified Bcsscl function. 

An excellent reference for the next step can be found 
in [30]. In order to generate a distribution of with p(0) 
given by (C5), one first generates a uniform deviate x 
(independent uniformly distributed random number be- 
tween and 1) and makes use of the fundamental trans- 
formation law of probabilities, which simply tells us 



\p(0)d0\ = \dx\. 
So we need to solve 



(C7) 



(C8) 



The solution of this is x — F(0), where F(0) is the indef- 
inite integral of p(0). The desired transformation which 
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takes a uniform deviate into one distributed as p(8) is 
therefore 



9{x) = F-\x) 



(C9) 



where F~ x is the inverse function to F. This process is 
illustrated in Figure 19. 

Unfortunately, F (and F~ 1 ) can only be computed nu- 
merically. In order to implement the method we used 
look-up tables and interpolation. On systems where in- 
teger operations are much faster than floating point op- 
erations, things can be speeded up considerably by dis- 
cretizing the 9ij (for instance one can take the integers 
to 524288 to correspond to to 2n) and then storing all 
possible values of the sinusoidal functions that can occur 
(all 524288 values). This requires some storage capacity 
(about 64 Mbyte for our implemention) but this should 
not be onerous for any machine that one would consider 
doing such simulations on. One should note that some 
machines can compute trigonometric functions in only a 
few clock cycles and therefore it may be faster than a 
look-up call to memory. The resulting code took about 
twice as long per Monte Carlo step (MCS) to run as the 
simple "pick at random and then reject" method, but this 
loss is more than compensated for by the orders of mag- 
nitude improvement in correlation times. There is still 
considerable freedom in the order in which subsequent 
lattice sites are selected. Naively, one would think that, 
as long as all sites are visited on some pseudo regular 
basis, that the order is unimportant. While this is true 
in the sense that the order is unimportant for eventually 
reaching equilibrium, the order can have a huge impact 
on how fast you get there. The slowest (in the sense of 
long correlation times) method is to select sites at ran- 
dom. One can significantly reduce (by a factor of up to 
about L depending on temperature) correlation times by 
going through the lattice in typewriter fashion or a mix- 
ture of random and typewriter ordering. However, one 
must go through in different directions (alternate left- 
right-up-down with up-down-left-right etc.) in order for 
the correlation times to be isotropic (i.e. have the same 
correlation time for say Y measured in both the x and 
y direction). To ensure the accuracy of the implemen- 
tation, the code was tested against published results for 
the / = and / = 1/2 cases. 



2. Error Analysis 

Suppose we make N successive observations A ll ,n = 
1, • • • , N, of a quantity A in our simulation. If the dis- 
tribution of the fluctuations in A is Gaussian (this is not 
true for all the parameters measured), then the expecta- 
tion value of the square of the statistical error, which in 
this case is the variance, is 



<(<L4) 2 > 



1 N 



1 N 



N 2 



N N 



jj«A 2 ) - (Af) 



li \ (MAJ - (Af 



(A 2 ) {AY 



The autocorrelation function for A is defined as 
(A A,) (A) 2 



(a 2 ) (Ay 



(CIO) 



(Cll) 



where we associate the time with step fi. Note that 
0a(O) = 1 and 4>A(t) decays to zero as t — > oo. The 
autocorrelation time ta is defined as 



(C12) 



For an exponential relaxation, ta is the relaxation time, 
so that for times 3> ta, <pA(tfi) is very small. If 
JV > ta we can, therefore neglect the term involving 
n/N in Eq.(ClO) and one obtains 



<(<L4) 2 ) = —((A 2 ) 



(Af)(l + 2T A ). 



(C13) 



Thus, our N correlated measurements are equivalent to 
N / (1 + It a) independent measurements, something that 
must be taken into account when calculating errors. 

The concept of self-averaging (or lack of) is extremely 
important in correctly estimating errors from Monte 
Carlo simulations with disorder. Suppose we meas ure A 
and calculate it's statistical error using yj ((<5 A) 2 ) from 
Eq.(C13). If yJ({5A) 2 ) reduces to zero if L — > oo (and 
N/(l + 2ta) fixed) we say A exhibits self-averaging. If, 
on the other hand, y/ ((5 A) 2 } reaches an L-indcpcndent 
nonzero limit, we say A exhibits a lack of self-averaging. 
Random systems exhibit a lack of self averaging near the 
critical point [31]. In fact, the distribution of most quan- 
tities (over realizations of disorder) is not even Gaussian, 
making the use of yj ((SA) 2 ) as a measure of the statisti- 
cal error somewhat questionable. 

In calculating errors we make use of, among other 
things, the bootstrap resampling technique described in 
[32] and more compactly in [27]. From the set of data T> 
produced by our Monte Carlo simulation we calculate a 
set x of parameters such as the energy, order parame- 
ter, etc. Due to the random sampling, T>q is not a unique 
realization of the true parameters x true . With different 
initial conditions or other slight variations we could have 
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measured any of an infinite number of other realizations 
T>\, T>2, ■ ■ ■■ Although the set x is not the true one x true , 
we assume that the shape of the probability distribution 
Xi — x , is the same, or very nearly the same, as the shape 
of the probability distribution x^ — x true . This is not an 
assumption that x and x trMe are the same, it is just as- 
suming that the way in which random errors enter the 
simulation does not vary rapidly as a function of x trMe , 
so that x can serve as a reasonable surrogate. 

Suppose we have in some way obtained a set of equiv- 
alent realizations of our data. For each realization T>j we 
calculate the parameters Xj in the same way as we ob- 
tained x from T> . Each simulated measured parameter 
set yields a point Xj — x . If we simulate enough data 
sets we can map out the desired probability distribution 
for the parameter space. As mentioned above, this dis- 
tribution of parameters is not necessarily Gaussian so we 
require some means of defining what we mean by the sta- 
tistical error. We take the statistical error to be width of 
the confidence region that contains 68% of the data (i.e. 
the confidence region is defined by the interval Xq ± a 
where, given the set of realizations of the parameter x, 
68% of the Xj lie between xq — a and xq + a). In this 
way, if our distribution is Gaussian, our definition of the 
error is just the standard deviation, as one would want 
for compatibility with the standard case. 

It only remains to explain how we obtain "a set of 
equivalent realizations of our data". The bootstrap 
method [32] used the actual data set T> , with it's n — 
N I (1 + 2ta) "independent" data points, to generate any 
number of synthetic data sets T>j, with n data points. 
The procedure is simply to draw n data points at a time 
with replacement from the set T> . For the bond disor- 
dered systems this includes bootstrap resampling of the 
set of realizations of bond disorder, as well as bootstrap 
resampling of the data from an individual realization of 
disorder. The basic idea behind the bootstrap is that the 
actual data set, viewed as a probability distribution, is 
the best available estimator of the underlying probability 
distribution. 
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FIG. 1. (a) Fluxoid pattern for ground states of / = | and 
/ = | (Unit cells are marked by solid lines) . Domain wall flux- 
oid pattern for / = |: (b) shift-by-one wall, (c) shift-by-two 
wall, (d) herringbone wall, and (e) herringbone wall with a 
shift-by- two (a vortex is shown as a dark square). 




FIG. 2. Partition of the square lattice into staircases with 
the current flowing up or down the staircases. 




FIG. 3. Illustration of several possible bends and kinks in 
the different types of domain walls, (a) A 90 degree bend 
in a / = 1/3 shift wall showing change from shift-by-one to 
shift-by-two wall, (b) / = 1/3 shift-by-one wall branching into 
two herringbone walls, (c) Kink in a / = 1/3 shift- by-three 
wall accomplished by moving the vortex marked in plaid. 
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FIG. 4. Energy of a square domain of size L x L in a sys- 
tem with periodic boundary conditions of size 120 x 120 for 
/ = 2/5. The line is the fit -0.0268(25) +0.344797(68)L 
+0.301(1) In L -1.28(3)(L/120) 4 . The inset shows the resid- 
uals for a linear fit (stars) and the fit including the quadrapole 
corrections (diamonds) . 
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FIG. 5. Energy of a square domain of size 15 x 15 in a 
system with periodic boundary conditions of size x x x for 
/ = 2/5. The line is the fit 5.961081(7) -1.086(l)(15/a;) 4 . 
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FIG. 6. il(x,x') T J- 1 I(x,x'), the lattice "Green's" 
function, for f — 1/3 along a slice in the x-direction in a finite 
size system with periodic boundary conditions along the direc- 
tion of the slice. The line indicates a fit to A(lnx + ln(L — x)). 
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FIG. 7. / = 1/3 Binder's cumulant U vs T for L = 36 to 
L — 84 (smaller L shown as dotted lines). 
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FIG. 8. / = 1/3 x vs T for L = 36 to L = 84 and scal- 
ing collapse of this data (inset) where x = (T — T c )L l t" , 
V = xL~'< /l/ , v = 1, and 7 = \. 
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FIG. 9. Finite size scaling plots for / = |. (a) logarithmic 
derivative of M at T c vs L, (b) specific heat maximum (hol- 
low) and at T c (solid) vs L, (c) \ at T c vs L, and (d) M at 
T c vs L. 



15 




0.2 0.22 0.24 0.26 0.28 
k B T/J 




-10 12 

L 1/ "(T-T C ) 

FIG. 10. left:min(p fe+ , p fe _) versus ksT/J. Note that 
data from larger L are smaller: mm(p k+ , pk-) vanishes as 
L — > oo as indicated by the finite-size scaling plot (right) 
which shows a reasonable collapse for min(pfc+, pfc_) ~ L a ^ v 
with a/v = -0.20 ± 0.02. 
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FIG. 11. (a) Specific heat for L = 24 to L = 96. The 
dashed line indicates T c . Note that the shoulder which ap- 
pears for intermediate lattice sizes goes away for the two 
largest L. This makes the scaling of C not as good as for 
the other variables, (b) Scaling collapse of data shown in (a) . 
(c) Power law scaling found by Lee and Lee for smaller sys- 
tem sizes, applied to the data shown in (a). The logarithmic 
scaling shown in (b) gives a better collapse of the data. In 
particular the lower curve in (c) , corresponding to the scaled 
L = 96 data is separating from the pack. 
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FIG. 12. Scaling collapse of Y where x = (T — T c )L 1/u , 
y = YL^I" , v = 1, and /3 = |. Inset: raw data (solid and 
dotted), f T (dashed), and a\T - T c f (dot-dashed). 
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FIG. 15. (a) Specific heat verses L and (b) susceptibility 
versus L 2 . Errors are comparable to the symbol sizes. 
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FIG. 13. Free energy as function of the negative of the en- 
ergy per site for / = 2/5 (5 = 0). A constant has been added 
to the curves in order to separate them. 
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FIG. 14. Free energy barrier vs system size for / = I and 
5 = (squares), 0.05 (circles) and 0.10 (triangles). S is the 
bond disorder strength. 
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FIG. 16. Finite size scaling plots for / = g,5 = 0.1: (a) 
logarithmic derivative of M vs L, (b) C/fcs vs L, (c) \ vs L, 
and (d) dM/dK and dY/dK vs L 
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FIG. 17. Solid on solid interface. Overhangs and bubbles 
are ignored in the SOS model and interface configurations can 
be described in terms of integer-valued height variables mea- 
sured from the straight, T — configuration of the interface. 




FIG. 18. Two solid on solid interfaces. The interfaces have 
a negative binding energy causing them to want to stick but 
they cannot cross. This "no crossing" condition results in an 
entropic repulsion which pushes the interfaces apart at high 
enough temperature. Zk is the separation of the interfaces at 
the fc'th step and A k is the number of steps the interfaces 
take in the same direction at the fc'th step. 
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FIG. 19. Transformation method for converting a uniform 
deviate x into a random deviate 6 distributed according to 
the function p(9). 



domain wall type 



energy per unit length 
f = 1/3 / = 2/5 



herringbone-0 
herringbone- 1 
shift-by- 1 
shift-by- 2 
shift-by-3 



0.05673742 J 
0.19503538 J 
0.11419998 J 
0.16666667 J 



0.08611726 J 

0.15889929 J 
0.16612232 J 
0.14764859 J 



TABLE I. Domain wall energies for stable domain wall 
structures (i.e. walls which produce a vortex pattern consis- 
tent with SH/Sdj — for every 9j). The n in herringbone-n 
denotes the associated shift where n — is the standard her- 
ringbone. 
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